program mkhis ! -- makes histgroms from list-mode data file.
  !
  ! DESCRIPTION:
  !
  ! There are two sisters of this: mkhisq and mkhisqp. They all
  ! produces the same results, i hope. You would prefer mkhisq if it's
  ! true. If it's not, you should believe the results of mkhis because
  ! it has the most redundant checking code for discarding invalid
  ! list-mode data.
  !
  !
  ! LICENSE AGREEMENT:
  !
  ! You agree that we don't agree to anything.
  !
  !
  ! DISCLAIMER OF WARRANTIES:
  !
  ! This software is provided to you "AS IS," and the authors and
  ! contributors disclaim any warranty or liability obligations to you
  ! of any kind, implied or express, including but not limited to the
  ! warranties of merchantability, fitness for a particular purpose,
  ! and non-infringement.
  ! 
  use kuty
  use mkhis_common
  implicit none
  interface
     subroutine decd(pe,pd,tad1,tad2)
       integer*8,intent(in),value::pe
       integer*4,intent(out)::pd
       integer*4,intent(out)::tad1(3)
       integer*4,intent(out)::tad2(3)
     end subroutine decd
  end interface
  integer*2,allocatable,target::hists(:,:,:,:,:)
  integer*2,pointer::histp(:,:,:,:)
  integer*2,pointer::histd(:,:,:,:)
  integer*2,pointer::hist(:,:,:,:)
  integer nok,noverflow,nflip,nunderflow,nlost
  integer*4,target::tad(3,2)
  integer*4 pd,tad_max(3,2),tad_min(3,2)
  integer*4,pointer::t1,a1,d1,t2,a2,d2
  integer,allocatable::off(:)

  !
  call parse_cmdline

  !
  allocate(hists(nx,ny,nx,ny,2),stat=istat)
  histd=>hists(:,:,:,:,1)
  histp=>hists(:,:,:,:,2)

  call mkhis_init

  nok=0
  noverflow=0
  nflip=0
  nunderflow=0
  noverflow=0
  nlost=0
  t1=>tad(1,1); a1=>tad(2,1); d1=>tad(3,1)
  t2=>tad(1,2); a2=>tad(2,2); d2=>tad(3,2)              
  
  !
  allocate(off(ndh))
  do id=1,ndh
     off(id)=mgap*(id-1)
  end do
  
  ifile=0
  file: do
     select case(get_next_file(ifile))
     case(-1) 
        exit
     case(0)
        cycle
     end select

     n=fs/4
     if(n>nbuf) then
        if(nbuf>0) deallocate(buf)
        allocate(buf(n),stat=istat)
        if(istat/=0) stop MEMFAIL
        nbuf=n
     end if

     b1=1
     read: do ie=1,n/fbs+1
        b2=min(b1+fbs,n)
        read(1,iostat=istat) buf(b1:b2)    
        if(istat/=0) stop READFAIL
        nb=nb+(b2-b1+1)
        b1=b2+1
     end do read
     close(1)

     lo_found=.false.

     decode: do ie=1,n
        if(lo_found) then
           flg=get_flg(buf(ie))
           if(flg==0) then
              
              ! pair found
              lo_found=.false.
              npair=npair+1
              call decd(loc(buf(ie-1)),pd,tad(1,1),tad(1,2))
              
              tad_max=max(tad,tad_max)
              tad_min=min(tad,tad_min)
              !
              if(t1>=nt.or.t2>=nt.or.a1>=na.or.a2>=na) then
                 write(*,*) "*** WARNING: index overflow detected."
                 write(*,*) "ie,t1,a1,t2,a2=",ie,t1,a1,t2,a2
                 noverflow=noverflow+1
                 cycle
              end if
              if(t1<0.or.t2<0.or.a1<0.or.a2<0) then
                 write(*,*) "*** WARNING: index underflow detected."
                 write(*,*) "t1,a1,t2,a2=",t1,a1,t2,a2
                 nunderflow=nunderflow+1
                 cycle
              end if
              !
              if(pd==1) then
                 hist=>histp
              else
                 hist=>histd
              end if
              tad=tad+1 ! <<<
              
              if(d1<=ndh) then
                 if(d2<=ndh) cycle
              else
                 if(d2>=ndh+1) cycle
                 nflip=nflip+1
                 write(*,*) "*** WARNING: index flip detected."
                 write(*,*) "d1,d2=",d1,d2
                 CYCLE ! <<< dont know what to do but for safty
!!$                 t1=>tad(1,2); a1=>tad(2,2); d1=>tad(3,2)
!!$                 t2=>tad(1,1); a2=>tad(2,1); d2=>tad(3,1)
              end if
              
              if(d1>ndh.or.d2>nd) then
                 write(*,*) "*** WARNING: index overflow detected."
                 write(*,*) "d1,d2=",d1,d2
                 noverflow=noverflow+1
              end if
              if(d1<=0.or.d2<=ndh) then
                 write(*,*) "*** WARNING: index underflow detected."
                 write(*,*) "d1,d2=",d1,d2
                 nunderflow=nunderflow+1
              end if
              !
              ix1=(d1-1)*na+(na-a1+1)+off(d1)
              iy1=t1
              ix2=(d2-(ndh+1))*na+a2+off(d2-ndh)       
              iy2=t2
              if(ix1>nx.or.ix2>nx.or.iy1>ny.or.iy2>ny) then
                 write(*,*) "*** WARNING: index overflow detected."
                 write(*,*) "t1,a1,d1,t2,a2,d2=",t1,a1,d1,t2,a2,d2
                 write(*,*) "ix1,iy1,ix2,iy2=",ix1,iy1,ix2,iy2
                 noverflow=noverflow+1
                 cycle
              end if
              if(ix1<=0.or.ix2<=0.or.iy1<=0.or.iy2<=0) then
                 write(*,*) "*** WARNING: index underflow detected."
                 write(*,*) "ix1,iy1,ix2,iy2=",ix1,iy1,ix2,iy2
                 nunderflow=nunderflow+1
                 cycle
              end if
              hist(ix1,iy1,ix2,iy2)=hist(ix1,iy1,ix2,iy2)+1
              if(pd/=1) ndelay=ndelay+1
              nok=nok+1
           else
              nlost=nlost+1
              if(flg==1) then
                 lo_found=.true.
              else
                 lo_found=.false.
              end if
           end if
        else        
           if(get_flg(buf(ie))==1) then
              lo_found=.true.
           else
              nlost=nlost+1
           end if
        end if
     end do decode
  end do file
  
  call print_stats(nf,nb*4,nok,ndelay,&
       npair,nflip,noverflow,nunderflow,nlost)
  
  if(nf>0) call output(nx*ny*nx*ny,loc(histp),loc(histd))
  
contains
  
  integer function get_flg(e)
    integer*4,intent(in)::e
    get_flg=ishft(e,-30)
  end function get_flg
  
end program mkhis
